Each research question draws its own challenges which are unique in themselves. But:

Data ➡️ Cleaning ➡️ Fitness evaluation ➡️ Analysis

The SBDI4R package enables the R community to directly access data and resources hosted by SBDI.

Please refer to the package documentation for details on how to install it. Once installed the SBDI4R package must be loaded for each new R session:

library(SBDI4R)

Various aspects of the SBDI4R package can be customized.

• caching

• e-mail

• download reason

The Swedish Electrofishing Registry - SERS (Department of Aquatic Resources, SLU Aqua).

2.8 M observations starting in the 1950’s.

fq_str <- pick_filter("resource") 

Follow the instructions. Your choices here would have been ‘in3’ ▶️ ‘dr10’ (data resource 10 = SERS). Your variable fq_str will now contain a string ‘data_resource_uid:dr10’.

y1 <- 2008
y2 <- 2012
fq_str <- c(fq_str, paste0("year:[", y1, " TO ", y2,"]"))
# Note the square brackets are hard limits

Both filter strings (for data resource and for time period) will be treated as AND factors. For references on how to use the filters see the SBDI APIs documentation.

Using the function occurrences() we can now query for the observations fulfilling our filter. If you haven’t specified your email and the download reason in the sbdi_config() before, you need to pass this here.

xf <- occurrences(fq = fq_str)

# Remove what is not a species
xf$data <- xf$data[xf$data$rank == "species",]

# Simply summarise all records by data source 
table(xf$data$dataResourceName)

## 
## SLU Aqua  Institute of Freshwater Research Swedish Electrofishing Registry - SERS 
##                                                                             93200

You can quickly plot all the observations as a PDF file with the function ocurrence_plot(), one page per species:

Note that the plot is saved to a .pdf file in the current working directory. You can find that with getwd().

occurrences_plot(xf, "output/obsPlot.pdf", 
                 grouped=FALSE, 
                 taxon_level="species", 
                 pch='.')

There are many other ways of producing spatial plots in R. The leaflet package provides a simple method of producing browser-based maps with panning, zooming, and background layers:

library(leaflet)
# drop any records with missing lat/lon values
xfl <- xf$data[!is.na(xf$data$longitude) | !is.na(xf$data$latitude),] 
marker_colour <- rep("#d95f02", nrow(xfl))
# blank map, with imagery background
leaflet() |> 
  addProviderTiles("Esri.WorldImagery") |>
  # add markers
  addCircleMarkers(xfl$longitude, xfl$latitude,  
                   radius = 1, 
                   fillOpacity =.5, 
                   opacity = 1,
                   col=marker_colour,
                   clusterOptions = markerClusterOptions())

A quick summary over the years reveals a drop in number of records over time.

hist(xf$data$year, 
     breaks = seq(y1, y2), 
     xlab = "Year")

In the same way we can summarise the number of observations for each species, by common or scientific name.

sppTab <- table(xf$data$commonName)
sppDF <- as.data.frame(sppTab)
colnames(sppDF)[1] <- "species"
head(sppDF)

##           species Freq
## 1                   61
## 2 Alpine bullhead 4615
## 3 American burbot 7081
## 4        Aral asp    6
## 5     Arctic char   46
## 6    aurora trout  856

sppTab <- table(xf$data$scientificName)
sppDF <- as.data.frame(sppTab)
colnames(sppDF)[1] <- "species"
head(sppDF)

##                                species Freq
## 1       Abramis brama (Linnaeus, 1758)   61
## 2   Alburnus alburnus (Linnaeus, 1758)  660
## 3   Anguilla anguilla (Linnaeus, 1758) 2140
## 4     Astacus astacus (Linnaeus, 1758)  618
## 5 Barbatula barbatula (Linnaeus, 1758)  620
## 6     Blicca bjoerkna (Linnaeus, 1758)   74

Perhaps, you want to send this table as a .CSV file to a colleague. Save the table:

write.csv(sppDF, "output/SERS_species_summary.csv")
# NOTE: again this will be saved on your working directory

Let’s now ask: How does the species richness vary across Sweden?

library(sf) # the function coordinates() and proj4string() are in sp
# load some shapes over Sweden's political borders
data("swe_wgs84", package="SBDI4R", envir=environment())
# a standard 50 km grid
data("Sweden_Grid_50km_Wgs84", package="SBDI4R", envir=environment())
grid <- Sweden_Grid_50km_Wgs84

obs <- st_as_sf(as.data.frame(xf$data),
                coords = c("longitude","latitude"),
                crs = st_crs(4326))

# overlay the occurrence data with the grid
ObsInGridListID <- st_intersects(grid, obs)
ObsInGridList <- lapply(ObsInGridListID, function(x) st_drop_geometry(obs[x,]))
wNonEmpty <- unname( which( unlist(lapply(ObsInGridList, nrow)) != 0) )
if(length(wNonEmpty)==0) message("Observations don't overlap any grid cell.")

The result ObsInGridList is a list object with a subset of the data for each grid cell. Now summarise occurrences within grid cells:

# apply a summary over the grid cells 
nCells <- length(ObsInGridList)
res <- data.frame("nObs"=as.numeric(rep(NA,nCells)),
                  "nYears"=as.numeric(rep(NA,nCells)),
                  "nSpp"=as.numeric(rep(NA,nCells)),
                  row.names = row.names(grid),
                  stringsAsFactors = FALSE)

cols2use <- c("scientificName", "year")
dataRes <- lapply(ObsInGridList[wNonEmpty], 
                  function(x){
                    x <- x[,cols2use]
                    colnames(x) <- c("scientificName", "year")
                    return(c("nObs" = length(x[,"scientificName"]),
                             "nYears" = length(unique(x[,"year"])),
                             "nSpp" = length(unique(x[,"scientificName"]))
                             )
                           )
                    }
                  )
dataRes <- as.data.frame(dplyr::bind_rows(dataRes, .id = "gridID"))
res[wNonEmpty,] <- dataRes[,-1]
resSf <- st_as_sf(data.frame(res, st_geometry(grid)))

And finally plot the grid summary as a map:

We may now ask whether species richness varies across latitude. So we go further by arranging the observations by latitude:

library(dplyr)
library(tidyr)
xgridded <- xf$data |>
    mutate(longitude = round(longitude * 4)/4, 
           latitude = round(latitude * 4)/4) |>
    group_by(longitude,latitude) |>
    ## subset to vars of interest
    select(longitude, latitude, species) |>
    ## take one row per cell per species (presence)
    distinct() |>
    ## calculate species richness
    mutate(richness=n()) |>
    ## convert to wide format (sites by species)
    mutate(present=1) |>
    do(tidyr::pivot_wider(data=.,  
                          names_from=species, 
                          values_from=present, 
                          values_fill=0)) |>
    ungroup()
## where a species was not present, it will have NA: convert these to 0
sppcols <- setdiff(names(xgridded),
                   c("longitude", "latitude", "richness"))
xgridded <- xgridded |> 
  mutate_at(sppcols, function(z) ifelse(is.na(z), 0, z))

And plot it accordingly:

library(ggplot2)

ggplot(xgridded, aes(latitude, richness)) + 
  labs(x = "Latitude (º)", 
       y = "Species richness") +
  lims(y = c(0,20)) +
  geom_point() + 
  theme_bw()

In this example we are interested in exploring opportunistically collected data from the Swedish citizen science species observation portal - Artportalen.

To begin with, we want be sure there is an unequivocal way to find the species within the order Odonata (dragonflies) and nothing else, so let’s search for “odonata”:

sx <- search_fulltext("odonata")
sx$data[, c("guid", "scientificName", "rank", "occurrenceCount")]

## [1] "https://species.biodiversitydata.se/ws/search.json?q=odonata&fq=idxtype%3ATAXON"

##       guid                          scientificName    rank occurrenceCount
## 1  9829523  Odonata associated gemycircularvirus 1 species               0
## 2 10072832  Odonata associated gemycircularvirus 2 species               0
## 3  8062407 Bdellodes odonata Wallace & Mahon, 1976 species               0
## 4      789                                 Odonata   order          207681
## 5  7367071    Ramalina fastigiata var. odonata Hue variety               0

We quickly see there that other taxonomic levels appear too, and also species that look suspiciously as not belonging to dragonflies. But there is only one order. Let’s refine the search. To know which search fields we can use to filter the search we use the function sbdi_fields(fields_type = "general"). The search field we are looking for is “order_s”.

sx <- search_fulltext(fq="order_s:Odonata", page_size = 10)
sx$data[, c("scientificName", "rank", "occurrenceCount")]

## [1] "https://species.biodiversitydata.se/ws/search.json?fq=order_s%3AOdonata&fq=idxtype%3ATAXON&pageSize=10"

##       guid                             scientificName  rank occurrenceCount
## 1  1429753                Gomphomacromia Brauer, 1864 genus               1
## 2  1426725               Austropetalia Tillyard, 1916 genus               0
## 3  4799335                 Sogjutella Pritykina, 1980 genus               0
## 4  4302686                    Neuragrion Karsch, 1891 genus               0
## 5  4799353              Xamenophlebia Pritykina, 1981 genus               0
## 6  1429769               Lauromacromia Geijskes, 1970 genus               0
## 7  1428195                     Sympetrum Newman, 1833 genus           27050
## 8  4798599 Corduliochlora Marinov & Seidenbusch, 2007 genus               0
## 9  1423625             Torrenticnemis Lieftinck, 1949 genus               0
## 10 1423468                  Cyanallagma Kennedy, 1920 genus               0

Now we can download the taxonomic data (note that the search is case-sensitive):

tx <- taxinfo_download("order_s:Odonata", 
                       fields = c("guid", "order_s","genus_s", "specificEpithet_s", 
                                  "scientificName",  "canonicalName_s", "rank"), 
                       verbose = FALSE)
tx <- tx[tx$rank == "species" & tx$genusS != "",] ## restrict to species and not hybrids

You can save the tx object as the complete species list for later use.

We start by searching for the data resource we are interested in using the function pick_filter(). This is an interactive query guiding you through the many resources available to filtering your query (data resources, spatial layers, and curated species lists).

# follow the instructions 
fq_str <- pick_filter("resource") 

Follow the instructions. Your choices here would have been “in3” –> “dr5”. Your variable fq_str will now contain a string “data_resource_uid:dr5”.

We only want to look at data from year 2000 to 2010:

y1 <- 2000
y2 <- 2010
fq_str <- c(fq_str, paste0("year:[", y1, " TO ", y2,"]"))
# Note the square brackets are hard limits

We also want to filter spatially for Southern Sweden (Götaland).

Vector spatial layers (eg. polygons) can be imported in a number of different ways. SBDI APIs take as search input polygons in the so-called WKT Well Known Text format. So the first step is to load a vector layer and transform it into a WKT string. You could instead use the data we provid in the SBDI4R package data("swe").

data("swe",package = "SBDI4R")
wGotaland <- swe$Counties$LnNamn %in% c("Blekinge", "Gotlands", "Hallands", 
                                        "Jönköpings", "Kalmar", "Kronobergs", 
                                        "Östergötlands", "Skåne", "Västra Götalands")
gotaland_c <- swe$Counties[wGotaland,]

There are details about this polygon that we need to take care before. The WKT string should not be too long to be accepted by the API service. Also, the polygon we just got is projected in the coordinate system SWEREF99 TM, and the API service only accepts coordinates in a geodesic coordinate system WGS84. Let’s construct the WKT string:

# transform the CRS
gotaland_c <- st_transform(gotaland_c,
                           crs = st_crs(4326))

# disolve the counties into one polygon
gotaland <- st_union(gotaland_c)

# create a convex hull of the polygon to simplify the geometry and 
# reduce the length of the WKT string
gotaland_ch <- st_convex_hull(gotaland)

# cast it as MULTIPOLYGON as this is what SBDIs API need
# NOTE: as of today, the SBDI APIs will only work properly if the polygon is 
# submitted as a MULTIPOLYGON
gotaland_ch <- st_cast(gotaland_ch, to = "MULTIPOLYGON")

# create WKT string
wkt <- st_as_text(gotaland_ch)

The WKT string then looks like this:

## [1] "MULTIPOLYGON (((13.33575 55.34003, 12.81633 55.38594, 11.25342 58.35786, 11.13161 58.90942, 11.13145 59.01184, 11.21142 59.0897, 11.31566 59.11651, 11.82032 59.23553, 11.94833 59.26237, 12.06197 59.27159, 12.23104 59.27357, 15.79383 59.03876, 15.84306 59.02498, 19.2889 57.99043, 19.3058 57.96888, 18.90037 57.44014, 18.86704 57.39753, 18.3725 57.00678, 18.30044 56.9528, 16.40805 56.20229, 14.19057 55.38557, 13.33575 55.34003)))"

Next, we download the observations using the command occurrences(), but be aware that the search fields may not be the same as those used to search for taxa. We therefore recommend using the function sbdi_fields("occurrence") to find out which search fields we can use to filter for occurrences. Here we see that the field we need this time is “order”.

xf <- SBDI4R::occurrences(taxon = "order:Odonata", 
                  fq = fq_str,
                  wkt = wkt,
                  extra = "collector")

We have now downloaded the data locally and depending on your configuration this will be cached on your computer. However, as the search and download could take long time, we recommend to save the data locally. appropriate

save(xf, file = "an_appropriate_name.rdata")
load(file = "an_appropriate_name.rdata")

Before we can use the observation records we need to know if the observation effort (sampling effort) has varied over time and in space. We can approximate observation effort from the data by defining field visits i.e. occasions at which an observer has sampled observations. We reconstruct field visits (that is, assign each observation a visitUID) using using the package BIRDS. Additionally we want the data to be summarized over a grid of 25 km (provided through the SBDI4R package). The following functions will perform many different summaries at the same time. Please refer to the BIRDS package documentation for more detail.

library(BIRDS)

OB <- organiseBirds(xf$data, sppCol = "species" , 
                    # We only want observations identified at the species level
                    taxonRankCol = "rank", taxonRank = "species", 
                    # the visits are defined by collector and named locality
                    idCols = c("locality", "collector"), 
                    timeCols = c("year", "month", "day"), 
                    xyCols =c("longitude","latitude") )

# We don't need the whole grid, just the piece that overlaps our searching polygon
wInt <- unlist(st_intersects(gotaland, Sweden_Grid_25km_Wgs84))
gotaland_grid25 <- Sweden_Grid_25km_Wgs84[wInt,]

SB <- summariseBirds(OB, grid = gotaland_grid25, spillOver = "unique")

Once summarised, we can see over space and for a few selected years how the number of observations is distributed:

maxC <- max(SB$spatial$nObs, na.rm = TRUE)
palBW <- leaflet::colorNumeric(c("white", "navyblue"), 
                               c(0, maxC), 
                               na.color = "transparent")
oldpar <- par()
par(mar = c(1,1,1,1), mfrow=c(1,3))
plot(SB$spatial$geometry, col=palBW(SB$spatial$nObs),
     border = "grey", main="All years") ## with palette
legend("bottomleft", inset = c(0,0.05),
       legend = round(seq(0, maxC, length.out = 5)),
       col = palBW(seq(0, maxC, length.out = 5)),
       title = "Number of \nobservations", pch = 15, bty="n")

## or export other combinations, e.g. one map per observed year
yearlySp <- exportBirds(SB, 
                        dimension = "spatial", 
                        timeRes = "yearly", 
                        variable = "nObs", 
                        method = "sum")

maxC <- max(yearlySp$'2005', na.rm = TRUE)
palBW <- leaflet::colorNumeric(c("white", "navyblue"), 
                               c(0, maxC), 
                               na.color = "transparent")

plot(yearlySp$geometry, col=palBW(yearlySp$'2005'), 
     border = "grey",main="2005")
legend("bottomleft", inset = c(0,0.05),
       legend = round(seq(0, maxC, length.out = 5)),
       col = palBW(seq(0, maxC, length.out = 5)),
       border = "grey",
       title = "Number of \nobservations", pch = 15, bty="n")

maxC <- max(yearlySp'2010', na.rm = TRUE)
palBW <- leaflet::colorNumeric(c("white", "navyblue"), 
                               c(0, maxC), 
                               na.color = "transparent")

plot(yearlySp$geometry, col=palBW(yearlySp$'2010'), 
     border = "grey",main="2010")
legend("bottomleft", inset = c(0,0.05),
       legend = round(seq(0, maxC, length.out = 5)),
       col = palBW(seq(0, maxC, length.out = 5)),
       border = "grey",
       title = "Number of \nobservations", pch = 15, bty="n")
par(oldpar)

We now want to use the number of field visits as the measure for sampling effort. :

library(cowplot)
library(ggplot2)
library(colorRamps)
library(gridExtra)

vis <- ggplot(data = SB$spatial, aes( fill = nVis)) +
  geom_sf() +
  ggtitle("Visits") +
  scale_fill_gradient(low = "#56B1F7",
                      high = "#132B43",
                      na.value = NA) +
  theme(plot.margin = margin(1, 1, 1, 1, "pt")) +
  theme_cowplot()

spp <- ggplot(data = SB$spatial ,aes( fill = nSpp))+
  geom_sf()+
  ggtitle("Number of species")+
  scale_fill_gradient(low = "#56B1F7",
                      high = "#132B43",
                      na.value = NA) +
  theme(plot.margin = margin(1, 1, 1, 1, "pt")) +
  theme_cowplot()

grid.arrange(vis, spp, ncol=2)

We see that SB contains an element called SB$temporal that contains a daily time series with time-specific rows when there is information. xts also supports day time, but dating below day resolution is not yet implemented in the BIRDS package.

sb.xts <- SB$temporal
head(sb.xts, 5)

##            nObs nVis nSpp
## 2000-03-24    1    1    1
## 2000-04-05    4    3    3
## 2000-04-06   11    6    3
## 2000-04-10    1    1    1
## 2000-04-12    3    3    1

Sub-setting is convenient in xts as you can do it with its dates and with a / for a range of dates.

sb.xts["2010-09-07"] #a specific day

##            nObs nVis nSpp
## 2010-09-07   19   10   12

sb.xts["2010-09-01/2010-09-15"] #for a period

##            nObs nVis nSpp
## 2010-09-01   46   19   14
## 2010-09-02   28   14   12
## 2010-09-03   23   10   10
## 2010-09-04   64   20   18
## 2010-09-05   74   27   12
## 2010-09-06   18    5   11
## 2010-09-07   19   10   12
## 2010-09-08   13    6    8
## 2010-09-09   32   12   14
## 2010-09-10    1    1    1
## 2010-09-11   16    9    8
## 2010-09-12   20   10    8
## 2010-09-13   14    5    9
## 2010-09-14    1    1    1
## 2010-09-15    3    3    2

sb.xts["2010-09"] #a specific month

##            nObs nVis nSpp
## 2010-09-01   46   19   14
## 2010-09-02   28   14   12
## 2010-09-03   23   10   10
## 2010-09-04   64   20   18
## 2010-09-05   74   27   12
## 2010-09-06   18    5   11
## 2010-09-07   19   10   12
## 2010-09-08   13    6    8
## 2010-09-09   32   12   14
## 2010-09-10    1    1    1
## 2010-09-11   16    9    8
## 2010-09-12   20   10    8
## 2010-09-13   14    5    9
## 2010-09-14    1    1    1
## 2010-09-15    3    3    2
## 2010-09-17    3    2    3
## 2010-09-18    9    5    5
## 2010-09-19   12    7    5
## 2010-09-21    3    2    3
## 2010-09-22    4    4    2
## 2010-09-23    3    3    2
## 2010-09-24   10    5    5
## 2010-09-25    7    4    6
## 2010-09-26    7    6    2
## 2010-09-28    2    2    2
## 2010-09-29    5    3    4
## 2010-09-30    2    2    2

The package xts has several tools for converting to different time periods. Here we use apply.monthly to obtain the total number of observations and visits per month. The plot command for an object of calss xts provides a many features. This makes it fairly easy to customize your plots. Read more in ?plot.xts.

library(xts)
obs.m <- apply.monthly(sb.xts$nObs, "sum", na.rm = TRUE)

# alternative for more summary options but slower
# obs.m <- exportBirds(SB,
#                      dimension = "temporal",
#                      timeRes = "monthly",
#                      variable = "nObs",
#                      method = "sum")

plot(obs.m, 
     col = "darkblue", 
     grid.ticks.on = "month", 
     major.ticks = "year", 
     grid.col = "lightgrey",  
     main = "Total number of daily observations and visits per month")

vis.m <- apply.monthly(sb.xts$nVis, "sum", na.rm = TRUE)
lines(vis.m, col = "orange", lwd=2, on=1)

We can now look at some particular species and ask whether those have changed in occurrence over time:

speciesSummary(SB)[,1:4]

##                       species nCells nObs nVis
## 1              Aeshna affinis      3   32   27
## 2             Aeshna caerulea      6   13   13
## 3               Aeshna cyanea    128  896  861
## 4              Aeshna grandis    156 1765 1735
## 5             Aeshna isoceles     25  172  169
## 6               Aeshna juncea    101  366  353
## 7                Aeshna mixta     81  677  646
## 8              Aeshna serrata     13   39   38
## 9           Aeshna subarctica     35  108   99
## 10             Aeshna viridis     13   40   35
## 11             Anax imperator     51  559  515
## 12            Anax parthenope      2    5    5
## 13        Brachytron pratense     93  505  492
## 14       Calopteryx splendens    102  682  628
## 15           Calopteryx virgo    124 1061 1007
## 16         Coenagrion armatum     23   74   67
## 17      Coenagrion hastulatum    118  941  908
## 18      Coenagrion johanssoni     15   75   70
## 19       Coenagrion lunulatum     42  111  104
## 20          Coenagrion puella    124 1415 1360
## 21      Coenagrion pulchellum    124 1439 1392
## 22     Cordulegaster boltonii     75  500  493
## 23             Cordulia aenea    122 1072 1053
## 24      Enallagma cyathigerum    149 1717 1630
## 25        Epitheca bimaculata     16   36   35
## 26           Erythromma najas     97  774  741
## 27       Erythromma viridulum     13  143  126
## 28      Gomphus vulgatissimus     48  160  155
## 29           Ischnura elegans    124 1351 1294
## 30           Ischnura pumilio     22  110   95
## 31               Lestes dryas     52  214  206
## 32              Lestes sponsa    147 1454 1385
## 33              Lestes virens     50  201  188
## 34     Leucorrhinia albifrons     47  185  180
## 35      Leucorrhinia caudalis     34  143  136
## 36         Leucorrhinia dubia     85  338  320
## 37    Leucorrhinia pectoralis     85  368  351
## 38     Leucorrhinia rubicunda     96  445  432
## 39         Libellula depressa    110  569  546
## 40            Libellula fulva      9  101   90
## 41   Libellula quadrimaculata    160 2096 2047
## 42        Nehalennia speciosa      3   34   33
## 43   Onychogomphus forcipatus     76  447  445
## 44      Orthetrum cancellatum    129 1079 1025
## 45     Orthetrum coerulescens     76  253  246
## 46       Platycnemis pennipes     92  545  529
## 47        Pyrrhosoma nymphula    127  974  945
## 48       Somatochlora arctica     28   45   42
## 49 Somatochlora flavomaculata     90  427  420
## 50     Somatochlora metallica    109  621  612
## 51             Sympecma fusca     50  269  264
## 52          Sympecma paedisca      2    5    5
## 53            Sympetrum danae    136  803  768
## 54        Sympetrum flaveolum     82  302  295
## 55     Sympetrum fonscolombii      2    2    2
## 56       Sympetrum sanguineum    137 1268 1224
## 57       Sympetrum striolatum     84  311  298
## 58         Sympetrum vulgatum    128 1072 1016

We pick two species and compare their trends in number of visits where the species where reported, relative to the total number of visits.

library(dplyr)
sppCountLq <- obsData(OB) |> 
            group_by(year,visitUID) |> 
            summarise("focalCount" = sum(scientificName == "Libellula quadrimaculata"),
                      "sppLength" = length(unique(scientificName)), 
                      .groups = "drop") |> 
            ungroup() |> 
            group_by(year) |> 
            summarise("focalCount" = sum(focalCount),
                      "nVis" = length(unique(visitUID)), 
                      .groups = NULL)

sppCountLq$relCount <- sppCountLq$focalCount/sppCountLq$nVis

sppCountSd <- obsData(OB) |> 
            group_by(year,visitUID) |> 
            summarise("focalCount" = sum(scientificName == "Sympetrum sanguineum"),
                      "sppLength" = length(unique(scientificName)), 
                      .groups = "drop") |> 
            ungroup() |> 
            group_by(year) |> 
            summarise("focalCount" = sum(focalCount),
                      "nVis" = length(unique(visitUID)), 
                      .groups = NULL)

sppCountSd$relCount <- sppCountSd$focalCount/sppCountSd$nVis

oldpar <- par(no.readonly = TRUE)
par(mar=c(4,4,1,1), las=1)
plot(sppCountLq$year, sppCountLq$relCount, type = "l", 
     lwd = 3, xlab = "Year", ylab = "Relative number of visits with observations", 
     ylim=c(0, max(sppCountLq$relCount)), xaxp=c(2000, 2010, 10))
lines(sppCountSd$year, sppCountSd$relCount, lwd=3, col="#78D2EB")
legend("bottomright", 
       legend=c("Libellula quadrimaculata","Sympetrum sanguineum"), 
       text.font = 3, col = c("black", "#78D2EB"), lwd = 3, bty = "n")

par(oldpar)